#include <unistd.h>     
#include <sys/types.h>  
#include <errno.h>      
#include <stdio.h>      
#include <stdlib.h>     
#include <pthread.h>    
#include <semaphore.h>  
#include <gmp.h>
#include <time.h>


void calc_a(void *ptr);
void calc_b(void *ptr);

sem_t mutex, mutex2,mutex3;
unsigned i = 0;
int stop = 0;

int main(int argc, char *argv[])
{

	pthread_t thread_a;
    pthread_t thread_b;
    mpf_t shared_variables[3];
	unsigned prec_bits = 35000000;
	time_t start, end;
	double dif = 0.0;
	sem_init(&mutex, 0, 0);
    sem_init(&mutex2, 0, 0);
    sem_init(&mutex3, 0, 0);  
	mpf_set_default_prec(prec_bits);
	mpf_init(shared_variables[0]);
	mpf_init(shared_variables[1]);
	mpf_init(shared_variables[2]);

    time(&start);
	pthread_create (&thread_a, NULL, (void *) &calc_a, (void *) shared_variables);
    pthread_create (&thread_b, NULL, (void *) &calc_b, (void *) shared_variables);
    pthread_join(thread_a,NULL);
    pthread_join(thread_b,NULL);
    time(&end);

	dif = difftime(end,start);
	printf("Time: %2lf\n",dif);
	printf("Iterations: %u\n",i);
	printf("PI = \n");
    mpf_out_str(stdout, 10, 10000001, shared_variables[0]);
    printf("\n");

    mpf_clear(shared_variables[0]);
    mpf_clear(shared_variables[1]);
    mpf_clear(shared_variables[2]);
    sem_destroy(&mutex);
    sem_destroy(&mutex2);
    sem_destroy(&mutex3);
	return 0;
}


void calc_a(void *ptr)
{
	mpf_t *y = ptr;
	mpf_t t;
	mpf_t temp;
	mpf_t aux;
	unsigned p;
	mpf_init(t);
	mpf_init(aux);
	mpf_init(temp);
	mpf_set_ui(*y,1); /* a = 1 */
	mpf_set_ui(t,4); /* t = 1/4 */
	mpf_ui_div(t,1,t);
	p = 1; /* p = 1 */
	sem_post(&mutex2);

	while(1)
	{
		sem_wait(&mutex);
		if(stop >= 1)
        {   
        	break;    
        }        
		mpf_set(aux,*(y+1));
		sem_post(&mutex2);
		mpf_add(aux,*y,aux);
		mpf_div_ui(aux,aux,2);

		mpf_sub(temp, *y, aux);
		mpf_pow_ui(temp, temp, 2);
		mpf_mul_ui(temp,temp, p);
		mpf_sub(t,t, temp);
		p <<= 1;
		mpf_set(*y,aux);
		sem_post(&mutex3);	
	}

		mpf_add(*y,*y,*(y+1));
		mpf_pow_ui(*y,*y,2);
		mpf_mul_ui(temp,t,4);
		mpf_div(*y,*y,temp);

		mpf_clear(t);
		mpf_clear(temp);
		mpf_clear(aux);

}

void calc_b(void *ptr)
{
	mpf_t *y = ptr;
	mpf_t aux;
	mpf_init(aux);
	mpf_sqrt_ui(*(y+1),2); /* b = 1/sqrt(2) */
	mpf_ui_div(*(y+1),1,*(y+1));

	while(1)
	{
		sem_wait(&mutex2);
		if(stop >= 1)
        {   
			mpf_clear(aux);
            sem_post(&mutex);
            pthread_exit(0);
        }
        mpf_set(aux,*y);
		sem_post(&mutex);
		mpf_mul(aux, aux, *(y+1));
		mpf_sqrt(aux,aux);
		sem_wait(&mutex3);
		mpf_set(*(y+1),aux);
		i++;

		if(i == 22)
        {
            stop++;
        }
	}
}

